rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
source("functions.R")
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggpubr)
library(fixest)
library(broom)
library(emmeans)
library(flextable)
library(officer)

theme_set(theme_sjplot())

# regression setting
lss_df <- readRDS("data/data_lss.RDS")
dict_df <- readRDS("data/data_dictionary.RDS")

colnames(dict_df) <- paste0("dict_", colnames(dict_df))
dat <- bind_cols(lss_df, dict_df)

events <- yaml::read_yaml("dictionary/events.yml")

# temporal dummies
dat <- dat |> 
  mutate(
    post_accords = if_else(date >= as.Date("2020-08-13"), 1, 0),
    us_embassy   = if_else(date >= as.Date("2018-05-14"), 1, 0),
    gaza_war2023 = if_else(date >= as.Date("2023-10-07"), 1, 0)
  )

# country dummies
join_countries <- c("Bahrain", "Morocco", "Sudan", "UAE")
dat <- dat |> 
  mutate(join = if_else(Country %in% join_countries, 1, 0))

gcc_countries <- c("Bahrain", "Kuwait", "Oman", "Qatar", "Saudi Arabia", "UAE")
dat <- dat |> 
  mutate(gcc = if_else(Country %in% gcc_countries, 1, 0))

axis_countries <- c("Iran", "Iraq", "Lebanon", "Syria", "Yemen")
dat <- dat |> 
  mutate(axis = if_else(Country %in% axis_countries, 1, 0))

peacetreaty_countries <- c("Egypt", "Jordan")
dat <- dat |> 
  mutate(peace_treaty = if_else(Country %in% peacetreaty_countries, 1, 0))

# issue dummies
dat <- dat |> 
  mutate(Stability = dict_stability,
         Economy = dict_economy,
         Security = dict_security,
         International_law = dict_international_law)




#############
stargazer(as.data.frame(dat), type = "text")
ggplot() +
  geom_histogram(aes(x = dat$lss), color = "white", bins = 30)+
  labs(x = "LSS", y = "Frequency")

#################################### all country
# OLS
reg00 <- lm(lss ~ Country, data = dat)
mtable(reg00,
       summary.stats=c("sigma","R-squared","F","p","N"))
plot_model(reg00, type = "pred", terms = c("Country"), axis.labels = "", 
           title = "",
           axis.title = c("Country", "Tone of reporting (LSS)"))

# fixed effect
reg01 <- feols(lss ~ 1 | Country + year, data = dat)
summary(reg01)

reg02 <- feols(lss ~ Country | year, data = dat)
summary(reg02)

# coef plot
coef_df <- broom::tidy(reg02)
head(coef_df)

coef_df <- tidy(reg02, conf.int = TRUE)
ggplot(coef_df, aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  labs(x = "Effect on LSS (relative to baseline country)", 
       y = "Country",
       title = "Country effects on LSS (controlling for year)") +
  theme_minimal()

plot_model(reg02, type = "est", sort.est = TRUE)

# margin plot
emm <- emmeans(reg02, ~ Country)
plot(emm)

emm_df <- as.data.frame(emm)
ggplot(emm_df, aes(x = reorder(Country, emmean), y = emmean)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  coord_flip() +
  labs(
    title = "Marginal means of LSS by Country (adjusted for year)",
    x = "Country",
    y = "Estimated marginal mean of LSS"
  ) +
  theme_minimal()



#################### country dummies
# OLS
reg10 <- lm(lss ~ join + gcc + axis, data = dat)
mtable(reg10,
       summary.stats=c("sigma","R-squared","F","p","N"))

# fixed effect
reg11 <- feols(lss ~ join + gcc + axis + peace_treaty | year, data = dat)
summary(reg11)

# coef plot
coef_reg11 <- tidy(reg11, conf.int = TRUE)
coef_reg11 <- coef_reg11 |> 
  filter(term %in% c("join", "gcc", "axis", "peace_treaty"))

ggplot(coef_reg11, aes(x = reorder(term, estimate), y = estimate)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  coord_flip() +
  labs(
    title = "Effects of Country Groups on LSS (adjusted for year)",
    x = "Country Group",
    y = "Estimated Effect on LSS"
  ) +
  theme_minimal()

plot_model(reg11, type = "est", sort.est = TRUE)

# margin plot
emm_df_reg11 <- emmeans(reg11, ~ join + gcc + axis + peace_treaty) |> 
  as.data.frame()

emm_df_reg11 <- emm_df_reg11 |> 
  mutate(term = case_when(join == 1 ~ "Join", 
                          gcc == 1  ~ "GCC", 
                          axis == 1 ~ "Axis",
                          peace_treaty == 1 ~ "Peace treaty")) |> 
  filter(!is.na(term))

emm_summary <- emm_df_reg11 %>%
  group_by(term) %>%
  summarise(
    emmean = mean(emmean),
    lower.CL = mean(lower.CL),
    upper.CL = mean(upper.CL)
  )

ggplot(emm_summary, aes(x = reorder(term, emmean), y = emmean)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  coord_flip() +
  labs(
    title = "Marginal Means of LSS by Country Group (adjusted for year)",
    x = "Country Group",
    y = "Estimated LSS"
  ) +
  theme_minimal()

###################### DiD of AA
## Estimate "Was the agreement effective overall?" (average treatment effect)
# all countries
reg30 <- feols(lss ~ post_accords * Country | Country + year, data = dat)
summary(reg30, cluster = "Country")

reg31 <- feols(lss ~ i(Country, post_accords) + 
                 us_embassy + gaza_war2023 | Country + year, data = dat)
etable(reg31)
iplot(reg31, main = "Impact of Abraham Accords by Country")

country_labels <- c(
  "Bahrain", "Egypt", "Iran", "Iraq", "Jordan", "Kuwait",
  "Lebanon", "Libya", "Morocco", "Oman", "Palestine",
  "Qatar", "Saudi Arabia" = "Saudi", "Sudan", "Syria",
  "UAE", "Yemen"
)
iplot(reg31, main = "Impact of Abraham Accords by Country", labels = country_labels)


# jointed
reg32 <- feols(lss ~ post_accords * join | Country + year, data = dat)
summary(reg32, cluster = "Country")

reg33 <- feols(lss ~ post_accords * join + us_embassy + gaza_war2023 | Country + year, 
               data = dat)
summary(reg33, cluster = "Country")

# gcc
reg34 <- feols(lss ~ post_accords * gcc + us_embassy + gaza_war2023 | Country + year, 
               data = dat)
summary(reg34, cluster = "Country")

# axis
reg35 <- feols(lss ~ post_accords * axis + us_embassy + gaza_war2023 | Country + year, 
               data = dat)
summary(reg35, cluster = "Country")

# axis
reg36 <- feols(lss ~ post_accords * peace_treaty + us_embassy + gaza_war2023 | Country + year, 
               data = dat)
summary(reg36, cluster = "Country")


########## DiD plot
# country effect
coef_df_reg31 <- data.frame(
  term = names(coef(reg31)),
  Estimate = coef(reg31),
  conf.low = confint(reg31, cluster = "Country")[,1],
  conf.high = confint(reg31, cluster = "Country")[,2]
)

coef_df_reg31 <- coef_df_reg31 |> 
  filter(grepl("^Country::.*:post_accords$", term)) |> 
  mutate(
    Country = gsub("^Country::(.*):post_accords$", "\\1", term)
  )


group_info <- data.frame(
  Country = c("Bahrain","Egypt","Iran","Iraq","Jordan","Kuwait",
              "Lebanon","Libya","Morocco","Oman","Palestine",
              "Qatar","Saudi Arabia","Sudan","Syria","UAE","Yemen"),
  group = c("Join","Peace treaty","Axis","Axis","Peace treaty","GCC",
            "Axis","Other","Join","GCC","Axis",
            "GCC","GCC","Join","Axis","Join","Axis")
)

coef_df_reg31 <- coef_df_reg31 |> 
  left_join(group_info, by = "Country")

ggplot(coef_df_reg31, aes(x = reorder(Country, Estimate), y = Estimate, color = group)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "Impact of Abraham Accords on LSS by Country",
    x = "Country",
    y = "Estimated DiD Effect",
    color = "Country Group"
  ) +
  theme_minimal()


# dummy effect
get_interaction_df <- function(reg, group_name) {
  coefs <- coef(reg)
  ci <- confint(reg, cluster = "Country")
  
  df <- data.frame(
    term = names(coefs),
    Estimate = as.numeric(coefs),
    conf.low = ci[,1],
    conf.high = ci[,2]
  ) %>%
    filter(grepl("post_accords:", term)) %>%
    mutate(
      Country = gsub(".*:(.*):post_accords$", "\\1", term),
      Group = group_name
    )
  
  return(df)
}

df_join <- get_interaction_df(reg33, "Join")
df_gcc  <- get_interaction_df(reg34, "GCC")
df_axis <- get_interaction_df(reg35, "Axis")
df_peace <- get_interaction_df(reg36, "Peace treaty")

df_all <- bind_rows(df_join, df_gcc, df_axis, df_peace) |> 
  mutate(term = case_when(
    grepl("join", term, ignore.case = TRUE) ~ "Join",
    grepl("gcc", term, ignore.case = TRUE)  ~ "GCC",
    grepl("axis", term, ignore.case = TRUE) ~ "Axis",
    grepl("peace", term, ignore.case = TRUE) ~ "Peace",
    TRUE ~ term
  ))

ggplot(df_all, aes(x = reorder(term, Estimate), y = Estimate, color = Group)) +
  geom_point(position = position_dodge(width = 0.7), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.2,
                position = position_dodge(width = 0.7)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "Estimated DiD Effects by Country Group",
    x = "Country",
    y = "Estimated DiD Effect",
    color = "Country Group"
  ) +
  theme_minimal() 



##############table
DiD_country <- etable(reg31, reg33, reg34, reg35, reg36,
                      cluster = "Country",
                      tex = FALSE, 
                      digits = 3,
                      signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10))
print(DiD_country)
write_html(DiD_country, file = "output/DiD_country.html")

# word
DiD_country_ft <- etable(
  reg31, reg33, reg34, reg35, reg36,
  cluster = "Country",
  tex = FALSE,
  digits = 3,
  signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10)
)

df_DiD_country <- as.data.frame(DiD_country_ft)

colnames(df_DiD_country) <- ifelse(colnames(df_DiD_country) == "", " ", colnames(df_DiD_country))
df_reg <- tibble::rownames_to_column(df_DiD_country, var = "Variable")

ft_DiD_country <- flextable(df_DiD_country) |> 
  theme_vanilla() |> 
  fontsize(size = 9, part = "all") |> 
  autofit() |> 
  set_table_properties(layout = "autofit", width = 1)

doc_DiD_country <- read_docx() |> 
  body_add_par("DiD Regression Results", style = "heading 1") |> 
  body_add_flextable(ft_DiD_country)

print(doc_DiD_country, target = "output/DiD_country.docx")


############### event study
## Estimate "How did the effect change over time?" (dynamic treatment effect)
# trend of non-/joined counties based on year 2019
avg_trend <- dat %>%
  group_by(year, join) %>%
  summarise(mean_lss = mean(lss, na.rm = TRUE), .groups = "drop")

ggplot(avg_trend, aes(x = year, y = mean_lss, color = factor(join))) +
  geom_line(size = 1) +
  geom_point() +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "blue") +
  annotate("text", x = 2020.2, y = max(avg_trend$mean_lss, na.rm = TRUE),
           label = "Abraham Accords\n2020-08-13", hjust = 0, vjust = -0.5,
           color = "blue", size = 4) +
  scale_color_manual(values = c("0" = "gray40", "1" = "red"),
                     labels = c("Non-join countries", "Join countries")) +
  labs(title = "Average reporting tone (LSS) over time",
       subtitle = "Comparison of join vs non-join countries",
       x = "Year", y = "Average LSS", color = "Group") +
  theme_minimal(base_size = 14)

# joint country
reg40 <- feols(lss ~ i(year, join, ref = 2019) | Country + year, data = dat)
etable(reg40)
iplot(reg40, main = "Event-study: Effect of Abraham Accords on reporting tone",
      xlab = "Year", ref.line = 2019, col = 2, ci = TRUE)

# not join country
reg41 <- feols(lss ~ i(year, I(1-join), ref = 2019) | Country + year, data = dat)
etable(reg41)

df_join <- tidy(reg40, conf.int = TRUE) |> 
  mutate(group = "Join")
df_nonjoin <- tidy(reg41, conf.int = TRUE) |> 
  mutate(group = "Non-Join")

df_all <- bind_rows(df_join, df_nonjoin)
df_all <- df_all %>%
  mutate(year_num = as.numeric(str_extract(term, "\\d{4}")))

ggplot(df_all, aes(x = year_num, y = estimate, color = group)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = df_all$year_num, labels = df_all$year_num) +
  labs(title = "Event Study: Effect of Abraham Accords",
       x = "Year",
       y = "Estimated LSS effect") +
  theme_minimal()


######## issue
reg50 <- feols(lss ~ post_accords * join + 
                 post_accords * Stability +
                 post_accords * Economy +
                 post_accords * Security +
                 post_accords * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg50, cluster = "Country")


# dummies
reg51 <- feols(lss ~ post_accords * join + 
                 post_accords * join * Stability +
                 post_accords * join * Economy +
                 post_accords * join * Security +
                 post_accords * join * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg51, cluster = "Country")

reg52 <- feols(lss ~ post_accords * gcc + 
                 post_accords * gcc * Stability +
                 post_accords * gcc * Economy +
                 post_accords * gcc * Security +
                 post_accords * gcc * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg52, cluster = "Country")

reg53 <- feols(lss ~ post_accords * axis + 
                 post_accords * axis * Stability +
                 post_accords * axis * Economy +
                 post_accords * axis * Security +
                 post_accords * axis * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg53, cluster = "Country")

reg55 <- feols(lss ~ post_accords * peace_treaty + 
                 post_accords * peace_treaty * Stability +
                 post_accords * peace_treaty * Economy +
                 post_accords * peace_treaty * Security +
                 post_accords * peace_treaty * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg55, cluster = "Country")


# all countries
reg54 <- feols(lss ~ post_accords * Country + 
                 post_accords * Country * Stability +
                 post_accords * Country * Economy +
                 post_accords * Country * Security +
                 post_accords * Country * International_law +
                 us_embassy + gaza_war2023 | Country + year, data = dat)
summary(reg54, cluster = "Country")

# table
DiD_country_issue <- etable(reg51, reg52, reg53, reg55,
                      cluster = "Country",
                      tex = FALSE, 
                      digits = 3,
                      signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10))
print(DiD_country_issue)
write_html(DiD_country_issue, file = "output/DiD_country_issue.html")

# word
DiD_country_issue <- etable(
  reg51, reg52, reg53, reg55,
  cluster = "Country",
  tex = FALSE, 
  digits = 3,
  signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10)
)

df_DiD_country_issue <- as.data.frame(DiD_country_issue)
colnames(df_DiD_country_issue) <- ifelse(
  colnames(df_DiD_country_issue) == "", " ", colnames(df_DiD_country_issue)
)

df_DiD_country_issue <- tibble::rownames_to_column(df_DiD_country_issue, var = "Variable")

ft_DiD_country_issue <- flextable(df_DiD_country_issue) |> 
  theme_vanilla() |> 
  fontsize(size = 9, part = "all") |> 
  autofit() |> 
  set_table_properties(layout = "autofit", width = 1)

doc_DiD_country_issue <- read_docx() |> 
  body_add_par("DiD Regression Results by Issue", style = "heading 1") |> 
  body_add_flextable(ft_DiD_country_issue)

print(doc_DiD_country_issue, target = "output/DiD_country_issue.docx")


###### plot issue
# reg50
coef_df_reg50 <- data.frame(
  term = names(coef(reg50)),
  Estimate = coef(reg50),
  conf.low = confint(reg50, cluster = "Country")[,1],
  conf.high = confint(reg50, cluster = "Country")[,2]
)

coef_df_reg50 <- coef_df_reg50 |>
  filter(grepl("post_accords", term)) |>
  mutate(
    Issue = case_when(
      grepl("Stability", term) ~ "Stability",
      grepl("Economy", term) ~ "Economy",
      grepl("Security", term) ~ "Security",
      grepl("International_law", term) ~ "International_law")) |> 
  filter(
    grepl("Stability|Economy|Security|International_law", term))

ggplot(coef_df_reg50, aes(x = Issue, y = Estimate)) +
  geom_point(size = 3, color = "steelblue") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "Effect of Abraham Accords by Issue",
    x = "Issue",
    y = "Estimated DiD Effect"
  ) +
  theme_minimal()

# reg51 to reg53, and reg55
coef_df <- function(reg, group_name){
  df <- data.frame(
    term = names(coef(reg)),
    Estimate = coef(reg),
    conf.low = confint(reg, cluster = "Country", level = 0.9)[,1],
    conf.high = confint(reg, cluster = "Country", level = 0.9)[,2]
  )
  df <- df |> 
    filter(grepl("post_accords:.*:Stability|post_accords:.*:Economy|post_accords:.*:Security|post_accords:.*:International_law", term)) |> 
    mutate(
      Issue = case_when(
        grepl("Stability", term) ~ "Stability",
        grepl("Economy", term) ~ "Economy",
        grepl("Security", term) ~ "Security",
        grepl("International_law", term) ~ "International_law"
      ),
      Group = group_name
    )
  return(df)
}

df_join <- coef_df(reg51, "Join")
df_gcc  <- coef_df(reg52, "GCC")
df_axis <- coef_df(reg53, "Axis")
df_peace <- coef_df(reg55, "Peace treaty")

df_all <- bind_rows(df_join, df_gcc, df_axis, df_peace)

ggplot(df_all, aes(x = Issue, y = Estimate, color = Group)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(width = 0.5), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "DiD Effect of Abraham Accords on LSS by Issue and Country Group",
    x = "Issue",
    y = "Estimated DiD Effect",
    color = "Country Group"
  ) +
  theme_minimal()


# by country: reg54
coef_df_country_issue <- function(reg){
  df <- data.frame(
    term = names(coef(reg)),
    Estimate = coef(reg),
    conf.low = confint(reg, cluster = "Country", level = 0.9)[,1],
    conf.high = confint(reg, cluster = "Country", level = 0.9)[,2]
  )
  df <- df |> 
    filter(grepl("post_accords:.*:Stability|post_accords:.*:Economy|post_accords:.*:Security|post_accords:.*:International_law", term)) |> 
    mutate(
      Issue = case_when(
        grepl("Stability", term) ~ "Stability",
        grepl("Economy", term) ~ "Economy",
        grepl("Security", term) ~ "Security",
        grepl("International_law", term) ~ "International_law"
      ),
      Country = gsub("^post_accords:Country(.*):.*$", "\\1", term)
    )
  return(df)
}

coef_df54 <- coef_df_country_issue(reg54)

ggplot(coef_df54, aes(x = Issue, y = Estimate, color = Issue)) +
  geom_point(position = position_dodge(width = 0.7), size = 2.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(width = 0.7), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  facet_wrap(~ Country, scales = "fixed") +
  labs(
    title = "DiD Effect of Abraham Accords on LSS by Country and Issue",
    x = "Issue",
    y = "Estimated DiD Effect",
    color = "Issue"
  ) +
  theme_minimal() +
  scale_x_discrete(labels = NULL)


